AD-A1Q3  078 

UNCLASSIFIED 


GEORGIA  INST  OF  TECH  ATLANTA  CENTER  FOR  THE  ADVANCEM— ETC  F/G  20/11 
NUMERICAL  ANALYSIS  OF  DYNAMIC  CRACK  PROPAGATION!  GENERATION  ANO— ETC(U» 
JUL  81  T  NISHlOK A  *  S  N  ATLURI  N00014-78-C-0636 

GIT-CACM-SNA-9  NL 


KD  A103078 


Office  of  Naval  Research 
Contract  N00014-78-C-0636  NR  064-610 
Technical  f(ep«rt  Vo.  11  v 
Report  No.  GIT-CACM-SNA  9  '/"V 


NUMERICAL  ANALYSIS  OF  DYNAMIC  CRACK  PROPAGATION:  GENERATION 
AND  PREDICTION  STUDIES 


BY 


Q~ 
O 
C 3 


s' 

ooj  w.ut;  is 
la  unlimited. 


Center  for  the  Advancement  of  Computational  Mechanics 
School  of  Civil  Engineering 
Georgia  Institute  of  Technology 
Atlanta,  Georgia  30332 

81  8  18  075/ 


NUMERICAL  AN LAYS IS  OF  DYNAMIC  CRACK  PROPAGATION:  GENERATION 


AND  PREDICTION  STUDIES 


T.  Nishioka*  and  S.N.  Atluri** 


Center  for  the  Advancement  of  Computational  Mechanics 
School  of  Civil  Engineering 
Georgia  Institute  of  Technology,  Atlanta,  Ga.  30332 


Abstract : 

Results  of  "generation"  (determination  of  dynamic  stress-intensity  fac¬ 
tor  variation  with  time,  for  a  specified  crack-propagation  history)  studies, 
as  well  as  "prediction"  (determination  of  crack-propagation  history  for  speci¬ 
fied  dynamic  fracture  toughness  versus  crack-velocity  relationships)  studies 
of  dynamic  crack  propagation  in  plane-stress/strain  situations  are  presented 
and  discussed  in  detail.  These  studies  were  conducted  by  using  a  transient 
finite  element  method  wherein  the  propagating  stress-singularities  near  the 
propagating  crack-tip  have  been  accounted  for.  Details  of  numerical  procedures 
for  both  the  generation  and  prediction  calculations  are  succinctly  described. 

In  both  the  generation  and  prediction  studies,  the  present  numerical  results 
are  compared  with  available  experimental  data.  It  is  found  that  the  important 
problem  of  dynamic  crack  propagation  prediction  can  be  accurately  handled  with 
the  present  procedures. 

Introduction: 

For  dynamic  crack  propagation  in  finite  elastic  bodies,  the  interaction 
with  the  crack-tip  of  stress  waves  reflected  from  the  boundaries  and/or  emanated 
by  the  other  moving  crack-tip  plays  an  important  role  in  determining  the  inten- 
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sit \.  of  the  dynamic  singular  stress-field  at  the  considered  crack-tip.  Be¬ 
cause  of  the  analytical  intractability  of  such  elasto-dynamic  crack-propagation 
problems,  computational  techniques  are  mandatory.  A  critical  appraisal  of 
several  such  computational  techniques  was  made  by  Kanninen  [1]  in  1978.  Most 
of  the  finite  element  techniques  reviewed  in  [1]  use  conventional  assumed  dis¬ 
placement  finite  elements  near  the  crack-tip  and  hence  do  not  account  for  the 

/ 

known  crack-tip  singularity.  Moverover  in  these  techniques,  crack-propagation 
was  simulated  by  the  well-known  "node  release"  technique,  which,  as  discussed 
in  [1],  may  not  be  sufficiently  accurate.  The  literature  on  dynamic  finite 
element  methods  for  simulation  of  fast  fracture,  since  the  appearance  of  [1], 
has  been  reviewed  in  [2,3,4]. 

In  Refs.  [2,3,4],  the  authors  have  presented  a  "translating-singular ity" 
finite  element  procedure  for  simulation  of  fast  crack  propagation  in  finite 
bodies.  In  this  procedure,  a  singular-element,  wherein  the  analytical  eigen 
functions  for  a  propagating'  crack  in  an  infinite  domain  were  used  as  basis 
functions  for  assumed  displacements,  was  used  near  the  crack-tip.  In  simu¬ 
lating  crack-propagation,  this  singular  element  was  translated  by  an  arbitrary 
amount  AE  in  each  time-increment  At  of  the  time-integration  scheme.  During 
this  translation,  the  crack-tip  retains  a  fixed  location  within  the  singular 
element;  however,  the  regular  isoparametric  elements  surrounding  the  moving  sin¬ 
gular  element  deform  appropriately.  It  was  shown  [2,3,4]  that  the  above  finite 
element  method,  which  was  based  [2,3]  on  an  energy-consistent  variational 
principle  for  bodies  with  changing  internal  boundaries,  leads  to  a  direct  eval¬ 
uation  of  dynamic  K-factors  for  propagating  cracks.  Attempts  at  simplifying 
the  above  procedure,  by  employing  alternatively,  a  singular  element  with  only 
the  well-known  Williams'  e i gen- funct  ions  for  a  stationary  crack  being  used  as 
element  basis-functions,  or  distorted  triangular  isoparametric  elements  (the  so- 
called  "quarter-point  elements"),  in  place  of  the  above  described  singular-ele- 


ment,  were  made  in  [5].  However,  all  the  examples  presented  in  [2-5]  fall  into  the 
category  of  "generation  studies"  in  the  sense  decribed  earlier.  Specifically, 
results  for  finite-domain  counterparts  of  the  well-known  analytical  problems 
for  infinite  domains,  solved  by  Broberg,  Freund,  Nillsson,  Thau  and  Lui,  Sih  et 
al  (as  referenced  in  [2,3]),  were  presented  in  [2-5],  to  indicate  the  effects 


of  finite  boundaries,  and  stress-wave  interactions,  on  dynamic  crack-tip  stress- 
intensity,  in  these  problems. 

In  the  present  paper,  which  emphasises  the  "inverse"  or  "prediction"  prob¬ 
lem,  namely  the  determination  of  crack-tip  propagation  history  in  a  plane 
stress/strain  problem  for  a  specified  dynamic  fracture-toughness  versus  crack- 
velocity  relation,  the  following  topics  are  discussed:  (i)  a  synopsis  of 
the  mathematical  formulation  for  analysis  of  the  "generation"  problem;  (ii)  des¬ 
cription  of  the  details  of  analysis  of  the  "prediction"  problem;  (iii)  detailed 
description  and  discussion  of  the  numerical  results  of  both  the  "generation" 
and  "prediction"  studies  of  wedge-loaded  rectangular  double  cantilever,  and 
tapered  double  cantilever,  beam  specimens  for  which  experimental  data  has  been 
reported  by  Kalthoff  et  al  [6,7]  and  independent  numerical  results  have  been 
reported  by  Kobayashi  et  al  [8],  and  Popelar  and  Gehlen  [9].  The  present 
paper  ends  with  some  conclusions  and  a  discussion  of  the  open  questions  in 
numerical  analysis  of  fast  crack  propagation  in  realistic  metallic  structures. 
Synopsis  of  the  Formulation  of  "Generation"  Problem: 

Consider  two  instants  of  time  t^  and  t^  “  t^  +  At.  Assuming,  without  loss 
of  generality,  that  the  crack  propagation  is  in  pure  mode  I,  let  the  crack  lengths 

at  t^  and  t^  be  and  T.  =  +  AT,  respectively.  Let  the  displacements,  strains, 

11  1  2  2 

and  stresses  at  t^  and  t2  be,  respectively,  (u^,  r  ,  and  o  ),  and  (u^,  C i j ’ 

2 

and  ).  The  variables  at  time  t^  are  presumed  known.  It  has  been  shown  [2,3] 
that  the  variational  principle  governing  the  dynamic  crack  propagation  between 


t^  and  t^  can  be  written  as: 
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In  the  above,  is  the  domain,  and  s  £  the  external  boundary  where  time-de¬ 
pendent  tractions  are  prescribed,  at  time  t are  the  prescribed  tractions 

+  -2 

at  time  t,  at  s  ,  (=s  „)  as  well  as  at  I,;  T.  are  the  prescribed  tractions  at 
1  ol  at  i 

time  t,  at  s  .  as  well  as  at  AE+;  (  )+  indicates  the  upper  half  of  the  crack 
Lai 

face,  which  only  is  considered  in  the  present  mode  I  problem.  It  is  seen  that 

are  the  cohesive  forces  holding  the  crack-faces  together  at  time  t^.  Thus, 

it  is  seen  that  the  integrand  (o^ v^)+( 6u^)+  in  the  last  term  of  the  r.h.s.  of 

Eq.  (1)  corresponds  to  the  term  of  energy-release  rate  due  to  dynamic  crack 

> 

propagation.  The  Eq.  (1)  may  thus  be  viewed  as  a  virtual  energy-balance  re¬ 
lation  for  dynamic  crack-propagation,  and  hence  the  present  numerical  method 

based  on  Eq.  (1)  is  inherently  energy-consistent. 

11  2  2  2 
In  Eq.  (1),  (Ui,o.  J  are  known,  while  , eij *  and  ui^  are  t*ie  variables. 

Now,  Eq.  (1)  is  used  to  develop  a  finite  element  approximation  at  time  t2*  Thus, 
the  domain  is  discretized  into  a  finite  number  of  elements,  with  a  domain 
immediately  surrounding  the  crack-tip  being  treated  as  the  so-called  "singular 
element",  and  the  domain  V^-V^  maPPec*  by  the  well-known,  8-noded,  isopara¬ 

metric  elements.  In  the  singular-element  V^,  the  basis  functions  for  assumed 
dispalcements  are  the  crack-velocity  dependent  eigen-function  solutions  to  the 
elasto-dynamic  problem  of  crack-propagation  in  an  infinite  domain,  as  discussed 
in  this  paper. 

Note  that  at  time  in  the  present  mode  l  problem,  the  crack  tip  is  loc¬ 
ated  at  x^  =  +  i\Y,  and  hence  the  s ingu .ar-element  is  centered  at  x^  =  +  AT.. 

In  developing  the  equations  for  the  finite  element  mesh  at  t£>  it  is  seen  from 
Eq.  (1)  that  the  variation  of  ujj  and  uj  must  be  known  in  the  finite  element  mesh 


at  t„.  However,  ,  u^»  and  were  solved  for,  in  the  finite  element  mesh 

at  t^.  In  the  mesh  at  t^  the  crack-tip  was  located  at  x^  =  E^  and  hence  the 
crack  element  was  centered  at  E^.  Thus  between  t^  and  t ^  (t^+At)  the  crack 
element  is  translated  by  an  amount  AE.  While  the  crack-element  is  trans¬ 
lated,  only  the  elements  surrounding  the  moving  crack-tip  are  distorted.  Thus 
the  finite  element  meshes  at  times  t^  and  t ^  differ  only  in  the  location  of  the 
crack-tip  (and  hence  the  crack-element)  and  the  shapes  of  the  immediately  sur¬ 
rounding  isoparametric  elements.  Thus,  the  known  data  at  and  u^  in  the 
mesh  at  t^  is  interpolated  easily  into  corresponding  data  in  the  mesh  at  t£. 
Further  details  of  the  above  translating-singularity-element  method  of  simu¬ 
lating  dynamic  crack  propagation  in  arbitrary  shaped  finite  bodies  can  be  found 
in  [2,3]. 

We  now  remark  briefly  on  the  basis  functions  for  assumed  displacements  used 
in  the  singular  element.  Let  xa(a*l,2)  be  fixed  rectangular  coordinates  in  the 
plane  of  the  present  2-dimensional  elastic  body,  with  the  crack-tip  moving  along 
the  x^  axis  and  is  no~mal  to  the  crack-axis.  We  introduce  a  coordinate  sys¬ 
tem  (.KyX-y)  which  remains  fixed  w.r.t.  the  propagating  crack-tip,  such  that 
;  =  x^-vt,  where  v  is,  without  loss  of  generality,  the  constant  speed  of  crack- 
propagation.  It  can  be  shown  [2,3]  that  the  elastodynamic  equations,  governing 


this  problem,  for  the  wave  potentials  (dilatational)  and  <P  (shear)  are: 


?  2  2  2  2 


2  2 


2  2  2 


[l-(v/c  )‘](j  ]/JO  +  O  ]>/3xp  =  -(2v/c"j)(3$/3t30  +  (l/c‘)  (3‘ V3tfc) 


and  a  similar  equation  for  w,  except  that  c,  in  Eq .  (2)  is  to  be  replaced  by  c  , 

d  s 

where  and  c^  are  the  lilatational  and  shear  wave  speeds  respectively.  The 
"steady-state"  eigen-function  solution  to  the  homogeneous  part  of  Eq.  (2),  namely, 
the  solution  which  appears  time-invariant  to  an  observer  moving  with  the  crack- 
tip,  and  satisfies  the  prescribed  traction  conditions  on  the  crack  face  (f<0, 
x  =t0)  can  be  derived  easily,  as  indicated  in  [2]  and  elsewhere.  We  use  these 


eigen  function  solutions  for  an  infinite  body,  as  basis  functions  for  assumed 


displacements  within  the  "crack-tip-singularity-element".  However  to  satisfy 


the  full  Eq.  (2),  the  undetermined  coefficients,  0^.  below,  in  the  eigen 
function  expansion  are  taken  to  be  functions  of  time.  Thus,  within  the  singular 
element. 


uaU,x2,t)  =  uaj  U,x2,v)Bj  (t)  [a=l, 2;  j  =  l,2..N] 


(3) 


where  u  .  are  the  above  described  eigen-functions,  and  g.  are  undetermined 
aj  :  J 

parameters,  which  are  to  be  determined  from  the  finite  element  equations  for 
the  cracked  body. 

As  seen  from  Eq.  (3),  .the  eigen  functions  u^  depend  on  the  crack-tip 
velocity.  In  the  present  numerical  approach,  the  crack-tip  velocity  is  assumed 
to  be  constant  within  each  time-increment  At,  say  v^  between  t^  and  t^+At,  and 
v2  between  t2  and  t2+At,  etc.  Thus,  between  t^  and  t^+At,  the  eigen-functions 
embedded  in  the  singularity-element  correspond  to  velocity  v^  and  those  between 
t2  and  t2+At  correspond  to  velocity  v2 .  Thus,  the  present  finite  element 
method  is  capable  of  handling  non-uniform-velocity  crack  propagation. 

The  total  velocities  and  accelerations  of  a  material  particle  in  the  sin¬ 
gular  element,  within  each  time  step,  correspoinding  to  Eq.  (3),  can  be  written 
as: 


U  . . 


VU 

aJ ,  ■  J 


w 


and 

ii  =»  u  .0.  -  2vu  .B.  +  v^u  .  _rP.  (3) 

»  '*J  J  *j,  J  "J .  ■  J 

where  (),-  =  3()/3'.t  and  (  )  implies  a  time  derivative. 

The  salient  features,  pertinent  to  the  studies  reported  in  this  paper,  of 
the  present  method,  the  mathematical  details  of  which  are  reported  elsewhere 
[2,3],  are  as  follows: 

(i)  The  eigen  functions  u  ^  (u=l,2)  lead  to  the  familiar  (l//r)  singu¬ 
larities  in  strains  and  stresses.  Thus  the  coefficient  B  (t)  is  directly 
related  (to  within  a  scalar  constant)  to  the  dynamic  stress  intensity  factor. 


Kx(t). 


b 


(ii)  The  compatibility  of  displacements,  velocities,  and  accelerations  of 

material  particles  at  the  boundary  of  surrounding  elements  with  those  of  the 

surrounding  (usual)  isoparametric  elements  is  satisfied  through  a  continuous 

least  squares  approach.  If  the  displacements,  velocities,  and  accelerations 

of  the  nodes  at  the  boundary  of  the  singular-element,  V  ,  are  q,  q,  and  d  re- 

s 

spectively,  the  above  least-squares  technique  leads  to  linear  algebraic  relations 
between  the  sets  (q,  q,  q)and  (8,  8,  and  8)  where  B  are  undetermined  parameters 
in  the  eigen-function  expansion,  Eq.  (3),  in  the  singular-element.  From  these 
equations  and  the  final  finite  element  equations  governing  the  nodal  displace¬ 
ments,  velocities,  and  accelerations  of  the  cracked  structure,  the  variables 
8,  6,  8  can  be  computed  directly.  Thus,  the  dynamic  stress-intensity  factor, 
as  well  as  its  first  two  time  derivatives,  are  computed  directly  in  the  present 
procedure. 

(iii)  The  "transient”  finite  element  equations  are  integrated  in  time,  using 
the  well-known  Newmark's  8-method  [2,3]. 

(iv)  Because  of  the  use  of  the  eigen  functions  in  a  moving  coordinate  sys¬ 
tem,  as  in  Eq.  (3),  in  the  singular-element,  there  is  the  presence  of  an  "ap¬ 
parent"  damping  matrix  for  the  singular  element.  Further,  for  the  same  reason, 
this  damping  matrix  as  well  as  the  stiffness  matrix  of  the  singular-element,  are 
unsymmetric.  However,  the  stiffness  and  mass  matrices  of  the  surrounding  iso¬ 
parametric  elements  are,  of  course,  symmetric.  Thus  the  final  finite  element 
equation  system  will  have  a  "small"  degree  of  unsymmetry.  This  equation  system 
is  solved,  in  the  present  studies,  using  a  simple  iterative  scheme. 

Details  of  Analysis  of  Prediction  Problem: 

The  problem  hero  is  to  predict  the  time  histories  of  crack-length  [E(t>], 
crack-velocity  [£(t)  v(t)],  and  possible  crack-arrest,  for  a  specified  re¬ 
lationship  of  dynamic  fracture  toughness  [K^]  versus  crack-velocity  [ v ] .  Until 
verv  recently,  it  was  presumed  that  the  relationship  versus  v  was  a  unique 


material  property.  Recently,  however,  this  presumption  was  brought  to  question 
as  discussed  in  1 10] ,  due  to  Che  apparent  geometry  and  load-rate  dependence  of 
the  dynamic  fracture  toughness.  A  slight  specimen-geometry  dependence  of  dy¬ 
namic  fracture  toughness  versus  crack-velocity  relationship  was  noted  in  the 
experimental  results  of  Xalthoff  et  al  [6,7]*  Kanninen  et  al[10J  also  found 
that  dynamically  initiated  (impact  loading)  dynamic  crack-propagation  and  quasi- 
statically  initiated  dynamic  crack  propagation,  apparently  are  characterized 
by  markedly  different  toughness  properties.  Ways  out  of  this  apparent  impasse 
that  have  been  suggested  include:  (i)  to  postulate  the  dependence  of  dynamic 
fracture  toughness  on  second  (acceleration)  and  higher-order  time  derivatives 
of  the  crack-length  and  (ii)  considerations  of  nonlinear  effects  near  the 
crack-tip,  such  as  plasticity. 

The  problems  considered  in  the  present  paper,  however,  may  be  argued  to 
fall  into  the  reals  of  linear  elasto-dynamics.  Experimental  specimens  for  which 
the  present  analysis  is  applied,  made  of  Araldite  Bused  by  Kalthoff  et  al  [7] 
may  be  considered  to  be  effectively  linear-elastic,  eventhough  secondary  effects 
due  to  rate-dependent  viscoelastic  properties  of  the  specimen  may  be  present. 

In  any  event,  the  present  numerical  results  and  their  comparison  with  the  ex¬ 
perimental  data  may  effectively  serve  to  check  the  reasonableness  of  this  ap¬ 
proximation.  Further,  the  presented  analysis  procedure  can  easily  be  extended 
to  account  for  any  postulated  dependence  of  dynamic  fracture  toughness  on  crack- 


tip  acceleration  and/or  other  higher  order  time  derivatives  of  crack  length, 


ie  when  =  K^O': ,  F. ,  £.  . .  )  . 

With  this  motivation,  we  present  some  details  of  anlavsis  of  the  prediction 


problem  when  the  fracture  toughness  relation  is  given  in  the  form  =  K^(£). 

Thus,  this  analysis  cannot,  inherently,  either  add  to  or  lessen  the  controversy 


surrounding  the  specimen  geometry  dependence  of  K 


Let  the  prediction  problem  be  .  .ms iderod  to  have  been  solved  upto  time  t^. 
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In  order  to  find  the  solution  at  the  crack-velocity  at  t^t  namely, 

V2'^t2^  raust  156  f°und-  To  this  end,  it  is  first  noted  that  the  dynamic  stress- 
intensity  factor  can  be  written  as: 

Kr  =  K^t.v)  (6) 

Since,  in  the  present  procedure,  the  velocity  of  crack-propagation  is  assumed 
to  be  constant  within  each  time-step,  an  approximate  procedure  to  predict  the 

r 

velocity  at  [t^+(At)/2]  will  be  sought.  Using  double  Taylor  series  expansion, 
it  is  seen  from  Eq.  (6)  that: 

l’Vl 

(7) 

where  is  the  predicted  value  of  at  t^  +  (At/2).  One  can,  upon  expanding 
terms,  write  Eq.  (7)  as: 

KiP  =  WV  +  (T}  *i  (ti,vi)  +I(~2)2  *1  (tl’vl)  +  R 

=  01(t:1)  +  (y)  ^(tj)  +|(y)2  jjjUj)  +  R  =  K*p  +  R  (8) 

where,  (*)  =  9()/3t,  and  R  is  "residue"  of  the  Taylor  expansion  indicated  in 
Eq.  (8).  Note  that  use  is  made  of  the  salient  feature  of  the  present  analysis 
procedure,  that  B^(t)  =  Kj(t),  6^  being  the  coefficient  of  the  first  eigen-func¬ 
tions  as  in  Eq.  (3). 


Vh 


v^+Av)'  = 


i  /  A  fc  x  9  ,  A  3  i 

}  —  +  Av— } 


KjCt.v) 


Since  during  dynamic  crack  propagation,  K^.  =  K^.^,  using  the  predicted  KJp 
of  Eq.  (8)  and  the  specified  versus  E(t)  relation,  the  crack  velocity  v(s£) 

at  the  time  [t^  +  (At/2)]  can  be  predicted.  If  the  arrest  dynamic-toughness  is 

3  it  it  di  r  r 

,  crack-arrest  is  predicted  if  K^p  .  Thus,  in  the  present  procedure, 

crack  arrest  is  predicted  as  a  terminal  event,  if  any,  in  the  propagation  analysis 


Using  the  above  predicted  crack-velocity  value,  the  finite  element  system 


of  equations  at  time  t^,  based  on  Eq .  (1),  are  constructed,  and,  from  these, 
the  actual  dynamic  stress-intensity  factor  K^(t2)  [^B^(t^)]  is  computed.  Thus, 
the  actual  K^.  at  t^  +  (At/2)  is  computed,  as; 

W  ~  d/2)  |KI(t1)  +  KT(t2)  1  (9) 
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The  correlation  between  the  predicted  K  of  Eq.  (8)  and  the  actual  K^  of  Eq. 
(9)  can  be  seen  to  depend  on  the  residue,  "R" ,  of  Eq.  (8).  To  ensure  this  cor¬ 
relation,  a  further  approximation  is  introduced  in  the  present  work  that  the 
residue  R  at  t^  +  (At/2)  can  be  approximated  by  its  known  value  at  [t^-(At/2)], 
in  the  generic  sensed  Thus,  in  the  present  procedure,  the  generic  algorithm 

used  to  find  K  at  t,  +  (At/2)  can  be  written  as: 

Ip  1 

Vt+  r>  ■  W  +  <¥)2  6<C1>+  I<f>2  vv 


-  IV'fT'-Wf" 


(10) 


In  all  the  presently  reported  computations,  when  Eq.  (8)  with  R=0  was  used,  a 
maximum  error  of  the  order  of  3%  between  and  was  noted.  However,  Eq.  (10) 
was  used,  this  maximum  error  reduced  to  the  order  of  0.5%. 

We  now  discuss  the  "generation"  and  "propagation"  calculations  performed 
on  rectangular  as  well  as  tapered  double  cantilever  beam  specimens  of  Araldite 
B  materials.  The  results  are  compared  with  the  corresponding  experimental 
results  reported  by  Kalthoff  et  al  [6,7],  and  pertinent  conclusions  are  drawn. 
Generation  Calculations: 

To  demonstrate  the  "generation"  type  calculations,  we  first  treat  a  wedge- 
loaded  rectangular  double  cantilever  beam  specimen  (WL-RDCB) ,  the  crack-propagation 
histories  and  dynamic  stress-intensity  factor  histories  in  which  were  directly 
measured  by  Kalthoff  et  al  [6].  The  relevant  geometric  data  of  the  WL-RDCB  speci¬ 
men  are  indicated  in  Fig.  1  which  also  shows  the  finite  element  model  wherein 

the  moving-singularity-element  is  shown  hatched,  at  the  beginning  of  crack  propa- 

2 

gation.  The  material  constants  used  in  the  present  analysis  are:  E=3380  MN/m 
and  Poisson's  ratio,  v=0.33.  In  the  experiments  of  [6],  several  test  specimens. 


1)  As  can  be  expected  from  Eq .  (1),  inherent  numerical  errors  (usually,  very 
small)  in  the  present  formulation  are  oscillatory  in  nature  [5],  Thus  R(-H-At/2) 
is  approximated  by  -R(t-At/w). 
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wherein  cracks  were  initiated  from  blunted  notches  with  crack-propagation 


initiation  stress-intensity  factors  K_  larger  than  the  fracture  toughness  K  , 

Iq  Ic 

were  studied. 

Note  that  the  actual  loading  mechanism  in  the  experiment  is  closer,  in  num¬ 
erical  simulation,  to  loading  the  finite  element  model  at  point  A  in  Fig.  1,  with 
the  material  to  the  left  hand  side  of  line  BA  in  Fig.  1  also  considered  to  be 


participating  in  the  motion.  In  the  first  attempt  at  the  analysis,  however, 
the  loading  was  modeled  to  act  at  point  B  in  Fig.  1  instead,  and  the  material 
to  the  left  of  line  AB  was  not  modeled.  In  the  remainder  of  the  paper,  the 
numerical  model  wherein  load  was  applied  at  point  A  of  Fig.  1  and  the  material 
to  the  left  of  line  AB  (Fig.  1)  was  also  modeled,  is  often  referred  to  as  the 
"actual  loading  condition",  and  the  other  one  as  the  "simplified  loading  con¬ 
dition",  respectively. 

In  their  report,  Kalthoff  et  al  [6]  identify  the  RDCB  specimens  with 
3/2  •  3/2 

values  2.32  MN/ra  and  1.33  MN/m  ,  respectively,  as  specimens  No.  4  and  17. 
For  convenience,  the  same  identification  is  used  in  the  presently  reported 
numerical  simulation. 

As  noted  earlier,  the  "generation”  calculation  used  as  input,  the  experi¬ 
mentally  measured  crack  length  (and  hence  crack-velocity)  history.  The  output 
of  the  calculation  is  the  directly  computed  dynamic  stress-intensity  factor 
at  the  tip  of  the  propagating  crack  for  various  time  instants. 

Fig.  2  shows  the  considered  crack  velocity  and  length  history  for  RDCB 
specimen  4  as  reported  in  [6].  Fig.  2  also  shows  the  presently  computed  dynam¬ 
ic  K.j  as  a  function  of  time,  along  with  comparison  experimental  results  of  [6], 
and  numerical  results  of  Kobayashi  (8).  The  present  calculation  for  was 
performed  in  3  alternate  ways:  (i)  direct  computation,  since  is  same  as 
the  undetermined  parameter  in  the  element  basis  functions  as  mentioned 
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earlier,  (ii)  from  a  crack-tip  integral  which  gives  directly  the  crack-opening 

energy,  and  using  the  crack-velocity  dependent  relation  between  K^.  and  the 

energy-release  rate,  and  (iii)  calculating  fracture  energy  from  a  global  energy 

balance  relation.  It  is  seen  that  all  the  3  values  agree  excellently,  thus 

pointing  to  the  inherent  consistency  of  the  present  numerical  procedure.  It 
should  be  pointed  out  that  the  results  in  Fig.  2  were  based  on  using  the 

forementioned  "simplified  loading  condition".  As  seen  from  Fig.  2,  the  present 

numerical  results,  as  well  as  those  of  Kobayashi  [8].  exhibit  a  pronounced 

peak  as  compared  to  the  experimental  results,  even  thoug  the  peak  occurs  much 

later  in  the  present  results  as  compared  to  those  of  Kobayashi  [8] . 

Fig.  3  shows  variation  of  different  energy  quantities:  input  energy 

(W) ;  kinetic  energy  (T) ;  strain  energy  (U) ;  and  fracture  energy  (F) ,  for  RDCB 

specimen  4,  when  the  "simplified  loading  condition"  is  used.  It  should  be 

noted  that  in  the  present  procedure,  each  of  the  quantities  W,  T,  U,  and  F  is 

calculated  separately  and  directly.  Thus,  the  fact  that  U+T+F  is  equal  to  W 

at  all  times  (no  other  energy  dissipation  mechanisms  are  accounted  for  here) 

is  an  inherent  check  on  the  accuracy  of  the  calculation.  That  this  is  so  can 

be  seen  from  Fig.  3. 

The  Fig.  4  demonstrates  the  effects  of  the  alternate  loading-conditions 

employed  in  the  finite  element  model  of  RDCB  specimen  4.  In  both  the  cases, 

3/2 

the  model  is  loaded  so  that  K_  =2.32  MN/m  .  For  this  value  of  K_  ,  the  de- 

Iq  Iq’ 

formation  profiles  of  the  crack  face  when  the  load  is  applied  at  points  A  and 
B,  respectively,  are  shown  in  Fig.  4.  It  is  seen  from  Fig.  4  that  for  the  same 
value  of  K  :  load  (and  dispacements)  at  points  A  and  B,  respectively,  are: 

970. 7N  (and  0.615  mm)  and  972. 3N  (and  .74  mm).  Thus  when  the  load  is  modelled 
to  act  at  B  (the  so-called  "simplified  loading  case")  there  is  more  apparent 
input  of  energy  to  the  specimen  than  when  the  load  is  modeled  to  act  at  A 
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(the  so-called  "actual  loading  case")  when  an  identical  crack-length  history 

as  in  Fig.  3  is  used,  but  with  the  "actual  loading  condition",  the  computed 

dynamic  k-factors  are  shown  in  Fig.  5.  Comparing  Figs.  2  and  5,  it  is  seen 

that  an  apparently  small  modification  in  the  load-condition  modeling  contributes 

to  a  substantial  difference  in  the  k-factor  variation.  It  is  seen  that  the 
results  in  Fig.  5,  for  the  "actual  loading  case"  agree  remarkably  well  with  the 

experimental  results  (considering  the  possible  rate-sensitive  behaviour  of 

Araldite  B  as  opoosed  to  the  present  linear  elastic  modeling),  and  the  peak  in 

the  present  K-results  is  much  smaller  than  that  in  Kobayashis'  [8]  results. 

The  variation  in  energies  W,  U,  T,  F  for  the  "actual  loading  case"  is  shown  in 

Fig.  6.  Comparing  Fig.  3  and  6  it  is  seen  that  W  in  the  "simplified  loading 

case"  is  higher  than  in  the  "actual";  T  is  higher  in  the  "simplified"  than  in 

the  "actual",  and  that  the  variations  of  U  and  F  are  qualitatively  similar  in 

both  the  "loading  cases". 

The  effect  of  the  two  loading  cases  for  the  RDCB  specimen  17  is  exhibited 

3/2 

in  Fig.  7.  It  is  seen  that  for  the  same  value  of  “1.33  MN/m  ,  the  load 
(and  displacement)  at  points  A  and  B  are,  respectively:  556. 5N  (and  0.35  mm) 
and  557. 7N  (0.425  mm).  Thus,  once  again,  the  apparent  input  energy  to  the 
specimen  is  larger  in  the  "simplified  loading  case"  than  in  the  "actual 
loading  case".  This  anamoly  in  modeling  will  have  consequences  in  the 
"propagation"  or  "application"  phase  calculations  in  the  RDCB  No.  17  specimen 
to  be  discussed  later. 

"PROPAGATION"  (OR  "APPLICATION")  CALCULATIONS 

We  now  present  calculations  aimed  at  predicting  crack-propagation 
history,  and  possible  arrest,  given  tiie  initial  loading  conditions  and  using 
the  hypothesis  that  there  is  a  given  material  toughness  data  in  the  form  of  a 
dynamic  f ract ure- toughness- vers us-crack-veloc ity-relat ion.  Experimentally 


evidence  [6,7]  that  there  is  the  possibility  of  a  slight  geometry  dependence 

of  this  toughness  property.  The  material  toughness  data  surmised  from  the 

experimental  findings  of  [6,7]  for  RDCB  specimens,  and  tapered  double  cantilever 

beam  specimens  (TDCB)  are  shown  in  Fig.  8.  In  the  present  calculations,  the 

RDCB  and  TDCB  toughness  data  are  used  in  the  prediction  of  crack-propagation 

histories  in  RDCB  and  TDCB  specimens,  resperctively .  Calculations  based  on 
using  RDCB  toughness  data  for  analysing  TDCB  specimens,  and  vice  versa,  are 

not  orted  here. 

The  results  of  the  "propagation"  or  "application"  type  calculations  for 
RDCB  specimen  4,  using  the  toughness  property  data  of  Fig.  8  and  "simplified” 
boundary  conditions,  under  plane  stress  conditions,  are  shown  in  Fig.  9.  It  is 
seen  that  the  present  predicted  crack  length  at  arrest  is  larger  than  in  the 
experiment,  even  though  the  stress-intensity  factor  variation  correlates  well 
with  the  experimental  result  for  most  of  the  crack-propagation  history.  The 
respective  results  with  the  "actual  boundary  conditions",  and  under  plane  stress 
conditions,  are  shown  in  Fig.  10.  It  can  be  seen  from  Fig.  10  that  the  presently 
calculated  length  history,  crack-velocity  history,  as  well  as  the  K-factor 
variation,  are  all  in  remarkably  good  agreement  with  the  experimental  results. 

It  is  noted  that  the  present  peak  value  in  the  K-factor  is  much  closer  to 
the  experimental  result,  than  that  in  the  solution  by  Kobayashi  [8].  To  com¬ 
pare  the  effects  of  plane  stress  versus  plane  strain  conditions,  a  "propagation" 
calculation  was  performed  on  the  RDCB  No.  4.  Specimen,  with  "simplified  boundary 
conditions",  and  the  results  are  shown  in  Fig.  11.  Comparing  Fig.  9  and  11, 
it  is  seen  noticeable  difference  can  be  found  between  the  stress-intensity  fac¬ 
tor  variation  between  the  two  cases  in  the  initial  phase  of  the  crack  propagation 
history,  and  the  final  crack-arrest  length  is  much  higher  in  the  plane  stress 
case  than  in  the  plane-strain  case. 


14 


The  energy  variations,  U,  T,  F  and  W  for  the  cases:  (i)  plane-stress, 
simplified  loading  case,  (ii)  plane-stress,  "actual  loading  case",  and  (iii) 
plane-strain,  simplified  loading  case,  are  shown  in  Figs.  12,  13,  and  14  res¬ 
pectively.  Comparing  Fig.  12  and  13  it  is  seen  that  the  ratio  of  maximum 
kinetic  energy  to  input  energy  in  the  simplified  loading  case  (0.278)  is  much 

larger  than  in  the  actual  loading  case  (0.233),  while  the  crack  arrest  length, 
comparing  Figs.  9,  and  10,  is  much  larger  in  the  "simplified  loading  case"  than 

in  the  "actual  loading"  case.  Likewise,  comparing  Figs.  12  and  14,  it  is  seen 

that  the  ratio  of  maximum  kinetic  energy  to  input  energy  in  the  plane-strain 

case  (0.266)  is  smaller  than  in  the  plane-stress  case  (0.278),  while  the  crack 

arrest  length  is  smaller  in  the  plane  strain  case  as  compared  to  the  plane-stress 

case. 

The  crack-surface  deformation  profile  for  the  propagating  crack  in  RDCB 
No.  4  specimen  are  shown  in  Fig.  15  for  various  instances  of  time.  Noting  the 
essentially  linear  shapes  of  these  profiles,  except  asymptotically  close  to  the 
crack-tip,  the  possibility  exists  to  devise  simple  method  to  find  the  stress- 
intensity  factors  from  the  crack-mouth  opening  displacment.  This  possibility 
is  successfully  explored  in  [11].  Figs. 16  through  21  show  the  contours  of 
principal-stress  difference  values  at  various  instances  of  time  in  the  movine-sin- 
gularity  element  of  the  RDCB  4  specimen  model,  in  the  plane  stress  case.  The  sequen¬ 
tial  pictures  demonstrate  graphically,  not  only  the  singular-stress-field  but 
the  total  stress  field,  and  its  magnification  near  the  propagating  cr3ck-tip. 

Since  in  the  present  finite  element  method,  the  effects  of  stress-wave  inter¬ 
actions  are  accurately  accounted  for,  and  the  total  stress  (singular  as 
well  as  nonsingular)  field  can  he  computed  accurately,  results  similar  to  these 
as  well  as  the  results  for  cimcumferent ial  stress  (not  shown  here)  can  be 
used  in  the  analysis  of  crack-branching.  Such  studies  will  be  presented  else- 
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where,  shortly. 

The  results  for  RDCB  17  specimen  for:  (i)  plane  stress,  simplified 

loading  case,  and  (ii)  plane  stress,  actual  loading  case,  are  shown  in  Figs. 

22  and  23  respectively.  Comparing  Figs.  22  and  23,  it  is  seen  that  the  higher 

crack  arrest  length  in  the  simplified  loading  case  than  in  the  actual  loading 

case  can  be  attributed  to  the  higher  apparent  input  energy  in  the  former  than 

in  the  latter  case,  as  seen  from  Fig.  7.  In  the  plane-stress,  actual  loading 
case,  the  calculation  was  continued  for  a  sufficient  time  after  crack  arrest 

(t_320  sec),  and  the  observed  oscillation  in  K-factor  is  shown  in  Fig.  24.  This 
oscillation  is  qualitatively  similar  to  that  recorded  in  the  experiments  [6], 

The  results  for  the  plane-strain,  simplified  loading  case,  are  shown  in  Fig.  25. 
In  comparing  Figs.  22  and  25,  comments  essentially  similar  to  those  made  in 

comparing  Figs.  9  and  11,  can  be  made.  The  energy  variations  in  RDCB  17  specimen 

for:  (i)  plane-stress,  simplified  loading  case;  (ii)  plane-stress  actual 

loading  case,  and  (iii)  plane-strain,  simplified  loading  case,  are  shown  in  Figs. 
26,  27,  and  28,  respectively.  Again,  in  comparing  Figs.  26,  27,  and  28,  comments 
essentially  similar  to  those  in  connection  with  the  comparison  of  Figs.  12, 

13,  and  14,  respectively,  can  be  made.  Thus  there  is  a  correlation  between  the 
ratio  of  the  maximum  kinetic  energy  to  input  energy,  and  the  crack  arrest  length. 
The  crack-surface  deformation  profiles  at  various  instants  of  time  in  RDCB  17 

specimen  shown  in  Fig.  29  are  similar  to  those  in  Fig.  15  for  specimen  4.  Re¬ 

sults  such  as  in  Figs.  15  and  29  form  the  basis  for  methods  of  obtaining  K  from 
crack-mouth-opening  displacments  discussed  by  the  authors  elsewhere  [11]. 

The  finite  element  model  for  the  tapered  double  cantilever  beam  (TDCB) 
specimen  is  shown  in  Fig.  30.  The  cross-hatched  element  shown  in  Fig.  30  is 
the  authors'  moving  singularity  element,  and  the  mesh  shown  in  Fig.  30  is  thus 
at  the  beginning  of  crack  propagation. 
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As  in  RDCB  specimen,  two  loading  cases  were  considered:  (i)  the  edge 

loading  case  wherein  load  is  supposed  to  act  at  point  B,  and  (ii)  the  actual 

pin  loading  case  wherein  the  load  is  modeled  to  act  at  point  A  in  Fig.  30. 

Plane  stress  conditions  are  invoked  in  both  the  loading  cases.  The  influence 

of  loading  position  is  demonstrated  in  Fig.  31.  In  all  the  cases  shown  in 

3/2 

Fig.  31,  the  model  is  loaded  so  that  K^*2.08  MN/ra  .  As  the  loading  point 
approaches  to  the  specimen  surface  while  keeping  the  x^-coordinate  constant, 
the  displacement  at  the  loading  point  becomes  larger  while  the  reaction  force 
is  almost  constant,  thus  the  input  energy  to  the  specimen  becomes  higher.  The 
input  energy  in  the  edge  loading  (loading  point  B)  is  also  shown  in  Fig.  31. 

It  is  seen  that  the  input  energy  in  the  edge  loading  case  is  much  higher  than 
that  in  the  actual  loading  case. 

The  computed  results  for  K^t),  £(t)  and  £(t)  for  both  the  loading  cases 
are  shown  in  Figs.  32  and  33  respectively.  In  the  edge  loading  case  as  shown 
in  Fig.  32,  after  240  ysec  the  stress  intensity  factor  becomes  almost  constant 
and  the  crack  propagation  with  a  relatively  slow  speed  (v  =  100  m/sec).  In  the 
actual  loading  case,  however,  as  shwon  in  Fig.  33  the  crack  was  arrested  earlier 
than  in  the  edge  loading  case.  The  value  variation  with  crack  length  for 
the  actual  as  well  as  edge  loading  case  is  shown  in  Fig.  34.  The  result  in 
the  actual  loading  case  shows  better  agreement  with  the  experimental  results 
obtained  by  Kalthoff  et  al  [7]. 

The  energy  variations  for  the  edge  loading  and  actual  loading  cases  are 

shown  in  Fig.  35  and  36  respectively.  Comparing  Figs.  35  and  36  it  is  seen 

that  the  ratio  of  maximum  kinetic  energy  to  input  energy  in  the  edge  loading 

case  (0.132)  is  much  larger  than  that  in  the  actual  loading  case  (0.093).  As 

also  observed  in  the  RDCB  specimen,  the  ratio  of  maximum  kinetic  energy  to  input 

energy  correlates  with  the  crack  arrest  length,  i.e.,  increasing  I  with 

arr 
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the  increasing  value  of  (max  T/W) .  Here  the  correlation  between  the  total  energy 
(U+T+F)  and  the  input  energy  W  in  the  TDCB  specimen  is  much  better  than  in  the 
RDCB  specimen  as  shown  earlier. 

The  crack  opening  profiles  in  the  edge  loading  case,  at  various  instants 
of  time,  are  shown  in  Fig.  37.  Because  of  the  loading  at  the  edge  of  the 

t 

specimen,  these  profiles  are  distinctly  nonlinear  as  compared  to  those  in  actual 
loading  case  (see  Figs.  15  and  29  in  the  RDCB  specimen). 

Finally,  Figs.  38  to  42  exhibit  sequentially,  the  contours  of  principal- 
stress  difference  at  various  instants  of  time  in  the  moving-singularity  element 
of  the  TDCB  specimen  with  the  edge  loading.  It  is  noted  that  the  size  of  the 
moving-singularity  element  (16x8)  mm  for  the  TDCB  specimen  while  it  is  (42x21)  mm 
for  the  RDCB  specimen.  Comparing  Figs.  16  to  21  on  the  one  hand,  and  Figs.  38 
to  42  on  the  other,  it  is  seen  that  the  effects  of  crack-propagation  and  stress 
wave  interactions  are  more  complex  in  the  TDCB  specimen  than  in  the  RDCB  specimen. 
Concluding  Remarks: 

The  results  presented  above  indicate  that  the  presently  developed  compu¬ 
tational  procedures  are  capable  of  accurately  predicting  dynamic  crack  propa¬ 
gation  and  arrest,  based  on  the  hypothesis  that  there  exist  a  "reasonable" 
geometry  independent  material  property  in  the  form  of  a  dynamic- fracture-  tough- 
ness-versus-crack  velocity.  The  results  also  demonstrate  the  importance  of 
modeling  the  loading  conditions  and  other  boundary  conditions  highly  accurately, 
in  an  elastodynamic  crack  propagation  problem. 

However,  other  questions  that  may  be  germane  to  the  subject  of  dynamic 
fracture  mechanics  itself,  such  as  the  load-rate  sensitivity  of  dynamic  frac¬ 
ture  toughness,  etc. ,  need  to  be  resolved  before  the  power  of  the  present  pro¬ 
cedures  can  be  fully  tested.  These  questions,  while  not  forming  the  subject  of 
the  present  paper,  have  been  attempted  to  be  discussed  by  the  authors  [12],  and 
others  { 10 1  elsewhere. 
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Figure  Captions: 

Fig.  1:  Finite  element  model  for  RDCB  specimen 

Fig.  2:  Variation  of  dynamic  stress  intensity  factors  in  RDCB  No.  4  with 

simplified  loading  (generation  phase) 

Fig.  3:  Energy  variations  in  RDCB  No.  4  with  simplified  loading  (generation 

phase) 

Fig.  4:  Crack  opening  profiles  in  RDCB  No.  4  with  different  loading  conditions 

Fig.  5:  Variation  of  dynamic  stress  intensity  factors  in  RDCB  No.  4  with  actual 

loading  (generation  phase) 

Fig.  6:  Energy  variation  in  RDCB  No.  4  with  actual  loading  (generation  phase) 

Fig.  7:  Crack  opening  profiles  in  RDCB  No.  17  with  different  loading  con¬ 

ditions 


Fig.  8:  Crack  velocity  versus  dynamic  fracture  toughness  relations  for 

Araldite  B  epoxy  (Kalthoff  et  al) 

Fig.  9:  Variation  of  dynamic  stress  intensity  factors  in  RDCB  No.  4  with 

simplified  loading  (application  phase,  plane  stress) 

Fig.  10:  Variation  of  Dynamic  stress  intensity  factors  in  RDCB  No.  4  with 
actual  loading  (application  phase,  plane  stress) 

Fig.  11:  Variation  of  dynamic  stress  intensity  factors  in  RDCB  No.  4  with 
simplified  loading  (application  phase,  plane  strain) 

Fig.  12:  Energy  variations  in  RDCB  No.  4  with  simplified  loading  (application 
phase,  plane  stress) 

Fig.  13:  Energy  variations  in  RDCB  No.  4  with  actual  loading  (application  phase, 
plane  stress) 

Fig.  14:  Energy  variations  in  RDCB  No.  4  with  simplified  loading  (application 
phase,  plane  strain) 


Fig.  15:  Variation  of  crack  opening  profiles  in  RDCB  No.  4  (actual  loading) 
Fig.  16:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t«0.0usec) 

Fig.  17:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t»100usec) 

Fig.  18:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t-200psec) 

Fig.  19:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t*300usec) 

Fig.  20:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t”4 OOpsec) 

Fig.  21:  Contours  of  principal-stress  difference  in  RDCB  No.  4  (t"528usec) 
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Variation  of  stress  intensity  factors  in  TDCB  specimen  with  edge  loading 
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loading 
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TDCB  specimen 

Energy  variations  in  TDCB  specimen  with  edge  loading 
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